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Abstract 

We propose a new method to determine the dipole (and quadrupole) 
component of a distribution of cosmic ray arrival directions, which can be 
applied when there is partial sky coverage and/or inhomogeneous expo- 
sure. In its simplest version it requires that the exposure only depends on 
the declination, but it can be easily extended to the case of a small am- 
plitude modulation in right ascension. The method essentially combines 
a x 2 minimization of the distribution in declination to obtain the multi- 
polar components along the North-South axis and a harmonic Rayleigh 
analysis for the components involving the right ascension direction. 

1 Introduction 

The search of a dipolar anisotropy in the arrival directions of cosmic rays (CRs) 
is of crucial importance to better understand the CR origin, the distribution 
of their sources and the way they propagate through the magnetized galactic 
medium. Several analyses have looked for a first harmonic in the right ascen- 
sion distribution of the arrival directions, and some positive findings have been 
reported. These include some giving few per thousand anisotropies at energies 
~ 10 15 eV, others giving percent level anisotropies at ~ 10 16 eV (which are 
however in contradiction with more recent data from Kascade which obtains 
upper limits below these findings) [1], and also a 4% anisotropy at 10 18 eV 
has been reported by the Agasa collaboration [2] . Due to the small size of the 
anisotropies observed it is clear that these measurements are extremely delicate 
since non-uniformities in the sky coverage, due to e.g. asymmetries of the ex- 
perimental setup or to weather effects not properly accounted for, could mimic 
the signal searched. 

On the theoretical side, realistic models of CR diffusion in the Galaxy, in- 
cluding the drift component which actually turns out to give the dominant 
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contribution to the CR macroscopic currents for energies above that of the knee 
and up to the ankle [3], predict anisotropics steadily increasing with energy, 
with a typical size of 10~ 3 at 10 16 eV and reaching the percent level at 10 18 eV. 
These predictions of course depend on the details of the regular and turbulent 
galactic magnetic field models adopted, on the CR composition assumed as well 
as on the distribution of cosmic ray sources, which are believed to be galactic 
below the ankle. 

The Auger experiment will gather the largest CR statistics at EeV energies 
and above, and can hence study accurately the large scale anisotropies in these 
energy regimes. This should allow to test Agasa's findings at EeV energies and 
also to study the ankle region, which is of particular interest because it is there 
that a transition from a galactic to an extragalactic component is believed to 
occur, and hence the behavior of the anisotropies at these energies could provide 
important clues to understand this transition. 

Besides the 2D harmonic analysis, which determines the phase in right ascen- 
sion (in addition to the Rayleigh amplitude), it is clearly desirable to reconstruct 
the CR dipolc in 3D, i.e. getting its actual amplitude and direction in the sky 
(both in right ascension and declination). A procedure to achieve this in the 
case of an experiment with full sky coverage was put forward by Sommers some 
time ago [4] , and was recently generalized by Aublin and Parizot [5] to the case 
of partial exposure relevant for present experiments. The basis of Sommers' 
method is to take each event direction n,, with an associated exposure u>i, and 
obtain the dipole through 



with N = J2i ^i" 1 {i = 1) W labels the different events observed). This ap- 
proach (and its generalization to partial sky coverage as well, hereafter referred 
to as the SAP method) has the virtue of being very easy to implement once the 
exposure is estimated, but has the following possible drawbacks: 

• Since each event is weighted by the inverse of its associated exposure, 
different regions in the sky do not have their real statistical weight in the 
determination of the dipole, as they have all been artificially uniformised 
by dividing by the exposure. Moreover, events in low exposure regions 
have very large weights and can hence give rise to enhanced fluctuations 
in the results. 

• The uncertainties in both the amplitude and the direction of the dipole are 
not obtained from the data but rather using Monte Carlo simulations with 
a known dipole anisotropy and comparing the input dipole parameters 
with the reconstructed values. 

• The method has no means of diagnosing whether the departure from 
isotropy has a dipolar shape, and for instance an intrinsic quadrupole or 




(1) 



2 



even a localized excess in a region of the sky will lead to a non-vanishing 
dipole amplitude. In these situations the uncertainties estimated from the 
MC simulations will clearly have little or no meaning. 



2 The x 2 +Rayleigh method 

We will here consider an alternative method to determine a dipole signature 
which tries to overcome some of the difficulties just mentioned and apply it to 
simulated data sets. 

Let d be the dipole's direction, and 7 = acos(n • d) the angle between the 
dipole and the event's arrival direction h. A dipolar distribution should give 
rise to a CR flux (incident on the Earth) of the form 

d$ 

— oc (l + acos7), (2) 

where a is the dipole's amplitude (i.e. D = ad). When this flux is observed 
by an experiment with non- uniform exposure (integrated over time) u>(h), the 
expected event rates should behave as 

diV 

— =Ww(n)(l + acos7) (3) 

with the normalization, which depends on a and d, fixed to reproduce the total 
number of events observed. Suppose that now we try to fit a dipole signal along 
a direction a", along which the expected distribution would be 

A AT rift 

7 =W d9' w(n)(l + acos 7 ), (4) 

dcos7' J 

where 9' is the azimuthal angle around the axis a", cos 7' = dl ■ h and 

cos 7 = d ■ n — cos (3 cos 7' + sin (3 sin 7' cos(#' — 9' d ), (5) 

where cos/3 = d ■ d! and 6' d is the azimuthal angle, measured around d' , of the 
dipole vector. 

It is then clear that if u> were uniform in the sky, the term proportional 
to cos(0' — 6' d ) in the above integral would vanish, leading to the behavior 
dA^/dcos7' oc (1 + a cos (3 cos 7'), and hence the dipole's amplitude inferred, 
barring statistical fluctuations, would be ctcos/3, which is just the dipole com- 
ponent along the d! axis. One may then envisage that the dipole's direction 
could be obtained as the one maximizing the reconstructed dipole's amplitude, 
but however, for non-uniform exposures the cos(#' — 9' d ) term does not average 
to zero in general, so that this procedure could lead to a biased result. 
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There is however a particularly relevant case where this bias is absent, which 
is when one considers d! = z, i.e. the NS equatorial axis, since for the case of 
uniform exposure in right ascension (with an arbitrary declination dependence), 
the cos(#' — 6' d ) term in the integral will indeed vanish 1 , hence leading to a 
behavior dA r /dcos7 / oc (1 + a z cos 7'), with a z = a sin 6 being the amplitude of 
the z component of the dipolc. This then suggests to use a % 2 method to fit a 
dipolar distribution along the z direction to get an unbiased estimate of a z . 

To fit a dipole signal along a generic direction d', one can just take n 7 bins 
of cos 7' (e.g. n 7 = 10 is a reasonable choice, and we checked that using a larger 
value does not improve significantly the results), and compute the number of 
events in each bin Nj and the corresponding expected values 

where ASlj is the solid angle for which cos 7' falls in the j-th bin, with dN/dfl 
given by Eq. (3). We can then write the \ 2 function as 



where the sum is clearly restricted to the bins 2 with Nj > 0. One can then 
determine the value a m for which x 2 ( a ) 1S minimized, and obtain Aa such that 
X 2 (a m + Aa) = x 2 (a m ) + 1- 

We hence apply this method along the z direction to obtain first a z . Con- 
cerning the dipole component orthogonal to the z axis, the natural strategy 
is to apply Rayleigh's method. In particular, it has been noticed in [5] that 
this method gives, at least in the case of partial sky coverage, a better deter- 
mination of the dipole's right ascension than the SAP method. One can use 
then Rayleigh's method to get the best estimate of the dipole's right ascension 
and first harmonic amplitude. This method 3 consists simply of computing the 
quantities (here a is the right ascension) 

2 N 2 N 

A=—^2cosai and ^sina; (8) 

i=\ i=l 



1 The other situation in which the result is unbiased is of course when d! = d, since in this 
case /3 = 0. 

2 In case one bin ends up with a very small number of events, it may just be convenient to 
choose a different number of bins for the fit. 

3 The Rayleigh method is based on the assumption that the experimental exposure is uni- 
form in right ascension. This may require eventually to select a subset of the data corre- 
sponding to sidereal days in which this assumption holds to a specified level of accuracy or, as 
discussed further below, to generalise the method so as to include the effects of right ascension 
non uniformities, if these can be properly modelled. 
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and then get the first harmonic amplitude r and phase "J/ through 



\J A 2 + B 2 and * = atan 



B 



(9) 



In addition, as shown in rcf. [5], there is a simple relation between r and the 
original dipolc components a z and a± = a cos S in the case in which the exposure 
is independent of a, which is 



r = 



c 3 a±_ 



where 



Cl 



-L 



ci + c 2 a 2 



d<5 uj (S) cos S 



(10) 



C2 



C3 



dS co(5) cos <5 sin <5 



d5 lu(5) cos 2 5. 



(11) 



One can then reconstruct completely the dipole's direction and amplitude from 
the previous expressions. Moreover, once the dipole's direction is obtained, it 
is useful to redetermine the dipole's amplitude a using the \ 2 minimization 
method for that fixed direction, and in this way one also evaluates the \ 2 /dof 
for this fit, what provides a check of the dipolar shape of the data distribution. 



3 Results 

As an example, we show the results of applying this method to data sets of 
3 x 10 4 simulated events with an intrinsic dipole of 5% amplitude (a = 0.05) 
pointing towards (5, a) = (— 45°, 0°), assuming an experiment at the Auger 
location (latitude —35.2°) with an ideal geometric exposure (uniform in right 
ascension and with a zenith angle dependence dN/d9 oc sin 9 cos 6, assuming a 
maximum zenith angle for the events analyzed of 60 degrees). 

Let us for instance describe in more detail the results of one particular sim- 
ulated data set. For this simulation the reconstructed dipole has an amplitude 
a = 0.059±0.013 pointing towards (5, a) = (-52°, -6°), which is - 10° degrees 
apart from the input direction of the simulations. Let us notice that the error 
Aa depends essentially only on the dipole's declination and the statistics at 
hand, scaling with the overall statistics as N -1 / 2 , as expected. We find indeed 
that, for the assumed detector's location coincident with the Auger latitude and 
for maximum zenith angles of 60°, the empirical expression 

Aa ~ J^- (1 + 0.6 sin 3 \6\) (12) 
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reproduces quite reasonably the results of many simulations performed with 
different input parameters. 

Fig. 1 shows the distribution of the events (solid line), the expectations of 
an isotropic sky (dashed line) and the best fit obtained through the method 
outlined (dotted line), as a function of cos 7 = d ■ h. 
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Figure 1: Histogram of events vs. cos 7 (solid), of isotropic expectations 
(dashed) and best dipole fit (dotted) for the simulation described in the text. 

The x 2 value of the best fit is here 13.1, and notice that the number of 
degrees of freedom is just n' 7 — 1, where n' 7 is the number of bins which are 
not empty (which can be less than n 7 due to the partial sky coverage, and 
moreover in this last case it will depend on the dipole's orientation). For this 
example n' 7 = 10, giving a value for x 2 /dof = 1.46, which is not unreasonable. 
A plot of the x 2 function for the best fit direction, as well as the fitting function 
X 2 ( a ) = X 2 ( a m) + ((a — a m )/Aa) 2 , are displayed in Fig. 2. 

Making this analysis for a set of 500 simulations like the one just described, 
we got the amplitude, right ascension and declination values (average and dis- 
persion) which are displayed in Table 1. It is apparent from this Table that both 
methods give average values quite close to the input ones, so that their biases 
are small, and the x 2 + Rayleigh method leads to dispersions 4 in the recovered 
values smaller by about 20% than SAP's method. This is probably related to 

4 Clearly part of the discrepancy between the input and reconstructed values is just due 
to statistical fluctuations in the simulated data, while another piece will be associated to the 
reconstruction method itself. 
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Figure 2: x 2 values for the direction maximizing the dipole's amplitude and 
fitting parabolic function. 



the possible large enhancement of the effects of statistical fluctuations in the 
low exposure regions near the boundaries of the observed sky, which affect the 
Aublin-Parizot generalization of Sommers' method for partial sky coverage. 
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SAP 


5.8 


1.7 


-44.6° 


19.1° 


-0.3° 


21.6° 


X 2 + Rayleigh 


5.4 


1.3 


-44.5° 


16.1° 


0.2° 


17.6° 



Table 1: Results from MC simulations (as described in the text) for the average 
values and dispersion in the reconstructed dipole's amplitude, declination and 
right ascension. The input values used were ai npu t = 5%, 5 = —45° and a = 0°. 

In the last row of Table 1 we display the values of the dipole amplitude 
obtained as \J a 2 z + a\, but for comparison the values of a which results from 
the x 2 minimization once the direction of the dipole has been fixed using a z , r 
and \&, lead to an average value (a m ) = 5.4 ± 1.4%, and the x 2 values of the 
corresponding fits are characterized by (x 2 /dof) = 0.88 ± 0.45. 
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Another advantage of the method proposed here is that it really exploits the 
dipolar shape of the signal searched. For instance, we performed 200 MC sim- 
ulations with 3 x 10 4 events each with a distribution consisting of an isotropic 
background plus a quadrupole, taken for dcfiniteness with symmetry of revolu- 
tion around an axis q, i.e. such that 



We adopted a quadrupole amplitude Q — 0.1 because this value leads to 
an inferred dipole a ~ 0.05, as in the examples discussed previously, and the 
orientation q adopted also points towards (5, a) = (— 45°, 0°). Applying the 
method introduced above to reconstruct a dipole leads to (a) = 6.3±1.1% (with 
(dm) = 5.3 ± 1.1%), but the key point is that the best fit results have associated 
values of {\ 2 /dof) — 3.5±1.3, which are quite poor, and hence this configuration 
cannot be mistaken with a dipolar one. Applying the SAP method to these 
simulations gives reconstructed dipole's amplitudes with (a) — 4.2 ± 1.5% and 
the associated uncertainty in a that would have been estimated for this case is 
just the same as given in Table 1, i.e. Aa ~ 1.7%. 

4 Incorporating right ascension modulation ef- 
fects 

A basic underlying assumption of the Rayleigh method is the uniformity of the 
exposure of the experiment with respect to right ascension. This however does 
not always hold, due e.g. to power or to communication failures, so that the 
strategy usually followed in order to apply this method is to first select a subset 
of the data corresponding to whole sidereal days in which the running conditions 
of the experiment under consideration were stable. This procedure can however 
significantly reduce the amount of data available. For instance, the Kascade 
Collaboration [6] was left with only ~ 20 % of their original data when this 
kind of selection was performed. 

We will here introduce a procedure which allows to take into account right 
ascension modulation effects, due to e.g. non-uniform running conditions of 
the experiment, making possible hence to use the whole statistics to obtain the 
amplitude and phase of the right ascension modulation of the incident cosmic 
ray flux. 

Suppose one has obtained, by properly taking into account the detector's 
down time [7], the right ascension dependent experimental exposure u(S, a). 
This allows then to introduce a declination only dependent average exposure 
through 




(13) 




(14) 
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It is now possible to generalize Rayleigh's method by introducing 

A=-^y cosai ,, v 11 . , 15 



and 



2 uj(Si) 

— > COS Oii — —z 

N 'Y u(Si,ai 
2 ^ . w(5i) 



where 



iV^i. (17) 

Notice that now the contribution of each event to A, B and to N is weighted 
by the factor u>(6i)/u)(6i,ai), which takes into account the right ascension non- 
uniformities, but being of order unity does not overweight the contribution of 
events in regions of relatively low exposure. 

By analogy with the standard approach, a modified Rayleigh amplitude and 
phase can now be introduced through 



* = atan-^ , f = V A 2 + B 2 . (18) 
A 

It is possible to generalize now eqs. (10) and (11) by identifying 

r 2 f d*.. , G>(6) 

A = — ail u>(d,a) — (d,a)cosa— - — - 

= ^"^o J dS u>(6) cos (5 J da cos a [l + a\\ sin S + a± cos 5 cos(a — ay)] 
2ir 

= — ®oc 3 a± cosa; d (19) 
N 

where we assumed that the cosmic ray flux is dipolar, i.e. following <J> = $o(l + 
ad, ■ n), and a« = a z while a± is the amplitude of the dipole vector in the xy 
plane, with ad being the right ascension of the dipole vector. 
Similarly one can find that 

27T 

B = -^-$oc 3 a± sin a d (20) 

and 

7V = 27r$o(ci+a||C 2 ) (21) 

with 



ci = J d5 uj(5) cos 8 
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C3 

One then finds that 



c 2 = J dS uj(5) cos S sin <5 

= J dS u){5) cos 2 5. (22) 



r = 



c 3 a±_ 



Ci + C 2 fl|| 



(23) 



while the Rayleigh phase turns out to be just the right ascension of the dipole 
a d . 



5 Including the quadrupole 

One of the main concerns in the reconstruction of a given multipolc from data of 
a partially covered sky, as for example the dipole reconstruction discussed at the 
beginning, is the possible mixing with the other multipolcs. For example, as we 
have shown, a non- vanishing quadrupole can lead to a dipolar signal when only 
a region of the sky is observed with an inhomogeneous exposure. A first control 
on the true dipolar character of the large scale anisotropy is given by the value 
of the x 2 /dof. A more careful analysis can be performed by reconstructing 
the next multipole, namely the quadrupole, and quantifying in this way their 
relative strength and the amount of leaking among them. 

We propose here an extension of the x 2 +Rayleigh method to reconstruct a 
dipolar + quadrupolar signal. This distribution would give rise to a CR flux 

d$ ^ A , _ „:„ i , _ _„r , _ _„j„:„_ Qzz . 2 



— = <&o I 1 + a z sin 5 + a x cos S cos a + a y cos S sin a H — — sin 5 



, Qxx 2 r 2 i Qyy 2 r • 2 r\ >- \ 

+ cos dcos a H — cos dsin a + Q xy cos o sin a cos a 

+Qxz cos S sin 8 cos a + Q yz cos 5 sin 5 sin a) , (24) 



where we have taken the z axis along the north pole direction, and the quadrupole 
tensor is symmetric and traceless. The observed flux, assuming for simplicity 
that the exposure only depends 5 on S, is given by 

do =A ^W (25) 

where the normalization, which now depends on the amplitude and orientation 
of the dipole and quadrupole, is fixed to reproduce the total number of events. 

The first step is again to fit the distribution of events in 5 integrated over 
a, taking advantage of the fact that the exposure is uniform in a. Now this 

5 The extension for u)(S, a) is straightforward along the lines proposed in the previous 
section. 
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distribution will be a function of not only a z , but also of Q zz , namely 

oc I 1 h a z smd H — sin o I . (26) 



d sin (5 

The values of a z and Qzz are obtained from a x 2 minimization similar to that 
performed in Eqs. (6,7) (notice that cos 7' = sin 6 for an axis along the z direc- 
tion). 

The next step is to use a generalized Rayleigh's method including higher 
order harmonics in order to obtain the rest of the dipole and quadrupole com- 
ponents. This is performed by computing the quantities 

2 N 2 N 

A = — cos a 4 , B = — sin a i . 

i=\ i=\ 

2 N 2 N 

C = — ^ cos ai sin ai , D = — sin 2 ai , 

2 iv 2 jv 

E = — ^ cos ai cos & , F =— ^ sin ai cos <5i . (27) 

i=l i=l 

These quantities can be easily related to the multipole coefficients in Eq. (24) 
using X)i /(&><*») = / dfl u)(S)d<fr/dft f(S,a). We thus obtain 

A = —j^— { a x cos S + Q-cz cos 8 sin 5) 
S = —j^— ( a y cos <5 + Qyz cos <5 sin 5) 

7T$ , 



C = -^-j-Q X y cos 2 (5 

£> = 27 ^° + a z sin <5 + Q Z z Q sin 2 5 - ^ cos 2 + -^Q yy cos 2 (5 

= — («x cos 2 (5 + Qjez cos 2 S sin 
F = —j^- {o-y cos 2 (5 + Q yz cos 2 (5 sin (28) 
Finally, the total number of events is related to <&o, a z and Q zz through 

N = 2n<P (l + a z shiS+^- (3 sin 2 S - fj^j . 
In the previous expression we have used the notation 



f(5) = ^ dS cos 5uj(S)f(S). 



11 



The remaining dipolar and quadrupolar coefficients are then obtained as 



a y = 



ICi 

K> 



— - (a cos 2 S sin S — E cos 5 sin 5) 
(^B cos 2 S sin 5 — F cos 5 sin S^j 



Qxz = cos S — A cos 2 S 

Q yz = -j^- (f cos S - B cos 2 6j 



Qxy 77- - CK, 



cos 2 5 



Q V y = = \ K,\D — 1 — a z sin d — (^4 sin 2 (5 — cos 2 

yy cos 2 s v s v 

Qxx — Qyy Qzzi (^^) 



with 



£1 - -4- = l + a z sin(5+%(3sin 2 (5-l), 

/C2 = cos (5 sin (5 cos 2 (5 — cos 2 S sin 5 cos 6. (30) 

As an example, we generated two datasets of 4 x 10 5 events each, one with 
a dipole and the other with a quadrupolc. For the first, we took a dipolc 
amplitude a = 0.1 pointing along (5, a) = (— 45°, 0°). Applying the method 
introduced above, we found from the \ 2 minimization 6 along the z-axis that 
a z = -0.085 ± 0.011 and Q zz = -0.018 ± 0.017. The generalized Rayleigh 
method then gave: a x = 0.075, a y = -0.005, Q xx = 0.007, Q yy = 0.011, 
Qxy = —0.004, Q xz — 0.002 and Q yz = —0.004, in perfect agreement with 
the input values, which correspond to (a x ,a y ,a z ) = (0.0707,0,-0.0707) and 
Qij = 0. Notice that applying to this same dataset the \ 2 + Rayleigh method 
introduced initially, i.e. without allowing for a quadrupole contribution, leads 
to a z = -0.075 ± 0.005, a ± = 0.072, with a = 0.104 ± 0.004, and pointing only 
2° away from the input direction. 

Regarding the dataset with an input quadrupolc, for which we took the 
expression given in Eq. (13), with Q = 0.1 and q also pointing towards (5, a) — 
(-45°, 0°), we got from the \ 2 minimization that a z — -0.014±0.010 and Q zz = 
0.045 ± 0.016. For the other components we obtained a x — 0.010, a y = —0.002, 
Qxx - 0.055, Q yy = -0.100, Q xy = 0.000, Q xz = -0.129 and Q yz = -0.011, 
in perfect agreement with the input values, which for the quadrupole assumed 
correspond to Q zz = Q/2, Q xx = Q/2, Q yy = -Q and Q xz = -3Q/2, with all 
other parameters vanishing. 

6 In this case, since we have to determine two parameters, it is convenient to take a larger 
number of bins of cos 7, e.g. n 7 = 20. 
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6 Summary 



In this work we have introduced a method to obtain large scale anisotropics 
in cosmic ray arrival directions in three dimensions. In its simplest version, 
the method allows to obtain the dipole vector describing the anisotropy of the 
incident flux, assuming that the exposure is uniform in right ascension, as ap- 
proximately holds for surface arrays. The z component a z is obtained from a x 2 
fit to the event distribution along the NS axis, and the perpendicular component 
is obtained with a 2D Rayleigh analysis in the orthogonal plane. We compared 
this method with the one recently introduced by Aublin and Parizot, showing 
that the dispersions in the values obtained are typically 20% smaller, and also 
the x 2 fit along the reconstructed dipolc's direction provides a measure of the 
agreement between the arrival direction distribution and a dipolar shape. 

We then generalized the method to the case in which a known modulation, 
e.g. induced by detector's down time, induces non-uniformities of the exposure 
with respect to right ascension. Finally, we extended the method to also include 
a quadrupole anisotropy, and showed with simulated datasets that this allows to 
recover the three components of the dipole and the 6 independent components 
of the quadrupole reliably. 

The methods here introduced should be applicable to detectors like Auger, 
and should help determine the large scale patterns of the cosmic ray distribution, 
which can provide crucial information about the origin of the cosmic rays and 
the way they propagate through the Galaxy. 
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